subroutine calculate_characteristics_normal(points,vertex2D,nVer2D,nInCat2D,&
    params,numbersEl,time,MN,maxR,maxZ)
    common /characteristics/ nPoints,nMn
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials
    integer elemsCount2D,typeTask,timeLevel,count,searchElement2D
    character*4 Ext
    real*8 eps2d,eps3d,eps3dNS,radiusRED,time,vertex2D,params,points,strength_normal,x,y,r,z,MN,&
    xM,yM,zM,xN,yN,zN,dx,dy,dz,Ex,maxR,maxZ
    logical file_exists
    dimension time(nt),vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d),points(3,nPoints),numbersEl(4,nodesCount2D),strength_normal(3),&
    MN(6,nMn),Ex(nt)

    if(nPoints.ne.0) then
        do j=1,nPoints
            r = dsqrt(points(1,j)**2+points(2,j)**2)
            if(r.eq.0.d0)r=1d-3
            z = points(3,j)
            nEl=searchElement2D(vertex2D,nVer2D,r,z,nodesCount2D,elemsCount2D)
            write(Ext,'(i4.4)') j
            open(500,file='.\results\SP_E_NX.'//Ext,mode='write')
            write(500,"(/,A)") '      t (mc)         ExN (mB)'
            open(501,file='.\results\SP_E_NY.'//Ext,mode='write')
            write(501,"(/,A)") '      t (mc)         EyN (mB)'
            open(502,file='.\results\SP_E_NZ.'//Ext,mode='write')
            write(502,"(/,A)") '      t (mc)         EyN (mB)'

            timeLevel=1
            write(Ext,'(i4.4)') timeLevel-1
            INQUIRE(FILE='.\calculations\Hfi\v2.'//Ext, EXIST=file_exists)

            do while(timeLevel.le.nt.and.file_exists)
                call calculateStrength(vertex2D,nVer2D,params,nInCat2D,nEl,&
                points(1,j),points(2,j),r,z,strength_normal,timeLevel,numbersEl,&
                maxR,maxZ,nodesCount2D,elemsCount2D,nMaterials2d)
                write(500,*)time(timeLevel)*1000.d0,strength_normal(1)*1000.d0
                write(501,*)time(timeLevel)*1000.d0,strength_normal(2)*1000.d0
                write(502,*)time(timeLevel)*1000.d0,strength_normal(3)*1000.d0
                timeLevel=timeLevel+1
                write(Ext,'(i4.4)') timeLevel-1
                INQUIRE(FILE='.\calculations\Hfi\v2.'//Ext, EXIST=file_exists)
            enddo

            close(500)
            close(501)
            close(502)
        enddo
    endif

    if(nMn.ne.0) then
        count = 20
        do j=1,nMn
            xM = MN(1,j)
            yM = MN(2,j)
            zM = MN(3,j)
            xN = MN(4,j)
            yN = MN(5,j)
            zN = MN(6,j)
            dx = (xN-xM)/count
            dy = (yN-yM)/count
            dz = (zN-zM)/count
            do timeLevel=10,nt
                Ex(timeLevel) = 0.d0
            enddo

            do k=1,count
                x = xM+dx/2.d0+(j-1)*dx
                y = yM
                nEl=searchElement2D(vertex2D,nVer2D,r,z,nodesCount2D,elemsCount2D)
                timeLevel=10

                write(Ext,'(i4.4)') timeLevel-1
                INQUIRE(FILE='.\calculations\Hfi\v2.'//Ext, EXIST=file_exists)

                do while(timeLevel.le.nt.and.file_exists)
                    r = dsqrt(x**2+y**2)
                    if(r.eq.0.d0)r=1d-3
                    z = 0.d0
                    call calculateStrength(vertex2D,nVer2D,params,nInCat2D,nEl,&
                    x,y,r,z,strength_normal,timeLevel,numbersEl,maxR,maxZ,&
                    nodesCount2D,elemsCount2D,nMaterials2d)
                    Ex(timeLevel) = Ex(timeLevel) + strength_normal(1)*dx
                    timeLevel=timeLevel+1
                    write(Ext,'(i4.4)') timeLevel-1
                    INQUIRE(FILE='.\calculations\Hfi\v2.'//Ext, EXIST=file_exists)
                enddo
            enddo

            write(Ext,'(i4.4)') j
            open(500,file='.\results\SP_dV_NX.'//Ext,mode='write')
            write(500,"(/,A)") '      t (mc)         dV(mV)'
            do timeLevel=10,nt
                write(500,*)time(timeLevel)*1000.d0,Ex(timeLevel)*1000.d0
            enddo
            close(500)
        enddo
    endif
end
!*************************************************************************
subroutine calculateStrength(vertex,nVer,params,nInCat2D,nEl,x,y,&
    r,z,strength_normal,timeLevel,numbersEl,maxR,maxZ,nodesCount2D,elemsCount2D,nMaterials2d)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials
    real*8 x,y,vertex,r,z,eNR,hz1,hz2,params,sigma1,sigma2,ValueInPoint2d,&
    strength_normal,Hfi1,Hfi2,Hfi3,Hfi4,Hfi,eNZ,hr1,hr2,maxR,maxZ,k1,k2,dhFI1,dhFI2
    integer elemsCount2D,timeLevel
    character*4 Ext
    dimension vertex(2,nodesCount2D),nVer(4,elemsCount2D),eNR(4),params(nMaterials2d),&
    nInCat2D(elemsCount2D),strength_normal(3),eNZ(4),numbersEl(4,nodesCount2D)

    write(Ext,'(i4.4)') timeLevel-1
    nFile = 503
    open(nFile,file='.\calculations\Hfi\v2.'//Ext,access='direct',recl=8,mode='read')
        
    
    do k=1,4
        !R
        i=nVer(k,nEl)
        if(vertex(1,i).eq.maxR)then
            n1=4
            n2=2
        else
            n1=3
            n2=1
        endif

        nEl1=numbersEl(n1,i)
        nEl2=numbersEl(n2,i)

        if(nEl1.eq.0.and.nEl2.ne.0)then
            sigma2=params(NInCat2D(nEl2))
            hz2=vertex(2,nVer(n1,nEl2))-vertex(2,nVer(n2,nEl2))
            read(nFile,rec=nVer(n1,nEl2))Hfi1
            read(nFile,rec=nVer(n2,nEl2))Hfi2
            eNR(k)=-(Hfi1-Hfi2)/sigma2/hz2
        elseif(nEl1.ne.0.and.nEl2.eq.0)then
            sigma1=params(NInCat2D(nEl1))
            hz1=vertex(2,nVer(n1,nEl1))-vertex(2,nVer(n2,nEl1))
            read(nFile,rec=nVer(n1,nEl1))Hfi1
            read(nFile,rec=nVer(n2,nEl1))Hfi2
            eNR(k)=-(Hfi1-Hfi2)/sigma1/hz1
        elseif(nEl1.ne.0.and.nEl2.ne.0)then
            sigma1=params(NInCat2D(nEl1))
            sigma2=params(NInCat2D(nEl2))
            hz1=vertex(2,nVer(n1,nEl1))-vertex(2,nVer(n2,nEl1))
            hz2=vertex(2,nVer(n1,nEl2))-vertex(2,nVer(n2,nEl2))
            k1=hz1/(hz1+hz2)
            k2=hz2/(hz1+hz2)
            read(nFile,rec=nVer(n1,nEl2))Hfi1
            read(nFile,rec=nVer(n2,nEl2))Hfi2
            read(nFile,rec=nVer(n1,nEl1))Hfi3
            read(nFile,rec=nVer(n2,nEl1))Hfi4
            eNR(k)=-(Hfi1-Hfi2)/sigma2/hz2*k1&
            -(Hfi3-Hfi4)/sigma1/hz1*k2
        endif


        !Z
        if(vertex(2,i).eq.maxZ)then
            n1=4
            n2=3
        else
            n1=2
            n2=1
        endif

        nEl1=numbersEl(n1,i)
        nEl2=numbersEl(n2,i)

        if(nEl1.eq.0.and.nEl2.ne.0)then
            sigma2=params(NInCat2D(nEl2))
            hr2=vertex(1,nVer(n1,nEl2))-vertex(1,nVer(n2,nEl2))
            read(nFile,rec=nVer(n1,nEl2))Hfi1
            read(nFile,rec=nVer(n2,nEl2))Hfi2
            read(nFile,rec=i)Hfi
            eNZ(k)=Hfi/vertex(1,i)/sigma2+(Hfi1-Hfi2)/sigma2/hr2
        elseif(nEl1.ne.0.and.nEl2.eq.0)then
            sigma1=params(NInCat2D(nEl1))
            hr1=vertex(1,nVer(n1,nEl1))-vertex(1,nVer(n2,nEl1))
            read(nFile,rec=nVer(n1,nEl1))Hfi1
            read(nFile,rec=nVer(n2,nEl1))Hfi2
            read(nFile,rec=i)Hfi
            eNZ(k)=Hfi/vertex(1,i)/sigma1+(Hfi1-Hfi2)/sigma1/hr1
        elseif(nEl1.ne.0.and.nEl2.ne.0)then
            sigma1=params(NInCat2D(nEl1))
            sigma2=params(NInCat2D(nEl2))
            hr1=vertex(1,nVer(n1,nEl1))-vertex(1,nVer(n2,nEl1))
            hr2=vertex(1,nVer(n1,nEl2))-vertex(1,nVer(n2,nEl2))
            k1=hr1/(hr1+hr2)
            k2=hr2/(hr1+hr2)
            read(nFile,rec=nVer(n1,nEl1))Hfi1
            read(nFile,rec=nVer(n2,nEl1))Hfi2
            read(nFile,rec=nVer(n1,nEl2))Hfi3
            read(nFile,rec=nVer(n2,nEl2))Hfi4
            dhFI1=(Hfi1-Hfi2)/hr1
            dhFI2=(Hfi3-Hfi4)/hr2
            eNZ(k)=(dhFI2+Hfi1/vertex(1,i))/sigma2*k1+(dhFI1+Hfi4/vertex(1,i))/sigma1*k2
        endif
    enddo
 
    strength_normal(1) = ValueInPoint2d(eNR,vertex,nVer,nEl,r,z,nodesCount2D,elemsCount2D)*x/r
    strength_normal(2) = ValueInPoint2d(eNR,vertex,nVer,nEl,r,z,nodesCount2D,elemsCount2D)*y/r
    strength_normal(3) = ValueInPoint2d(eNZ,vertex,nVer,nEl,r,z,nodesCount2D,elemsCount2D)
    close(nFile)
end
!*************************************************************************
real*8 function ValueInPoint2d(v,vertex,nVer,num,r,z,nodesCount2D,elemsCount2D)
    integer elemsCount2D
    real*8 v,vertex,psi1,psi2,psi3,psi4,r1,r2,z1,z2,r,z
    dimension v(4),vertex(2,nodesCount2D),nVer(4,elemsCount2D)
    r1 = vertex(1,nVer(1,num))
    r2 = vertex(1,nVer(2,num))
    z1 = vertex(2,nVer(1,num))
    z2 = vertex(2,nVer(3,num))
    ValueInPoint2d = psi1(r1,r2,z1,z2,r,z)*v(1)+&
                     psi2(r1,r2,z1,z2,r,z)*v(2)+&
                     psi3(r1,r2,z1,z2,r,z)*v(3)+&
                     psi4(r1,r2,z1,z2,r,z)*v(4)
end
!*************************************************************************
!*************************************************************************
!*************************************************************************